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Abstract 

A numerical simulation is carried out for static vortices in a long Josephson 
junction with an exponentially varying width. At specified values of the parameters 
the corresponding boundary-value problem admits more than one solution. Each 
solution (distribution of the magnetic flux in the junction) is associated to a Sturm- 
Liouville problem, the smallest eigenvalue of which can be used, in a first approxi- 
mation, to assess the stability of the vortex against relatively small spatiotemporal 
perturbations. The change in width of the junction leads to a renormalization of 
the magnetic flux in comparison with the case of a linear one-dimensional model. 
The influence of the model parameters on the stability of the states of the magnetic 
flux is investigated in detail, particularly that of the shape parameter. The critical 
curve of the junction is constructed from pieces of the critical curves for the dif- 
ferent magnetic flux distributions having the highest critical currents for the given 
magnetic field. 

1 Statement of the problem 

A Josephson junction model that takes into account the influence of the shape in the xy 
plane was considered in Refs. |U El E]- For a model junction with generatrices that vary 
by an exponential law E] (see Fig. |TJ) the basic equation for the phase <p(t, x) in the 
junction can be written in the form 

(f + a(f - if" + siny? + g(t,x) = , x e (0,1) , t>0. (1.1) 

The spatial coordinate x is normalized to the Josephson penetration depth Aj, and 
/ = (L + L s + L L )/\j. The time t is normalized to the plasma frequency (see, e.g., 
the monograph cited as Ref. [3]). The parameter a describes dissipative effects in the 
junction. An overdot denotes differentiation with respect to time t, and a prime denotes 
differentiation with respect to the spatial coordinate x. 
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The term g(t,x) = a[(p'(t,x) — h B ] in Eq. models the current caused by the 

exponential variation of the width of the junction. The dimensionless shape parameter 
a > is determined by the following expression (see Fig. ^ for notation): 




) 



where it is assumed that Wq > Wl- It is known that after suitable normalization the 



function ip(t, x) can be interpreted as the (dimensionless) magnetic flux along the junction. 
Then the value of h B is the (dimensionless) strength of the external (boundary) magnetic 
field along the y axis. For simplicity it is assumed below that the field h B is independent 
of time. 



Figure 1: Diagram of the Josephson junction of exponentially varying width. 
A detailed derivation of Eq. (jl.lj) is given in Ref . [2] . 

In the present paper we restrict consideration to a model junction with an in-line 
geometry. If the injection current density through the end L of the junction is denoted 
by 7 (which is assumed constant), the boundary conditions at x = and x = I, after the 
transition to the corresponding limits, take the form 



Because of the presence of a dissipative term oop the solution <p(t, x) of equation 
can lose energy and, consequently, for t — > oo it will approach the corresponding static 
distribution ip s (x) (to simplify the notation below we shall use the subscript s when 
necessary). This makes it necessary to investigate in detail the static distributions in the 
Josephson junction and their behavior upon variations of the model parameters. 



The boundary-value problem for the static distributions follows directly from relations 
dEU and 




(p'(t,0) = h B -h, (p'(t,l)=h B . 



(1.2) 



if" + sinv? + g(x) = , x G (0, 1) 
(0) - /i B + Z 7 = 0, 
(p'(l)-h B = 0, 



(1.3a) 
(1.3b) 
(1.3c) 



where the "geometric" current is given by the expression 



g(x) = a [(f'(x) - h B ) . 



(1.4) 
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We note that the solutions of the problem (|1.3|) depends not only on the spatial coordinate 
x but also on the set of parameters p = {I, a, hs, 7}, i.e., (p(x,p). When necessary below 
we shall indicate the dependence of the quantities on the parameters. 

Let ip s (x) be some static solution of equation i.e., a solution of the boundary- 

value problem (jl.3|) . With the goal of investigating the stability of the solution (p a {x), 
following Ref. [Hj, we consider a spatiotemporal fluctuation of the form 

<p{t,x) = <p a {x) +ee~ at / 2 Y, [e^"Vn(x) + e'^^x)] , (1.5) 

n 

which depends on the small parameter e. Substituting the expansion (jl.5|) into Eq. (jl.lj) 
and conditions (|1.2jl . to a first approximation in e we arrive at the Sturm-Liouville (SL) 
problem 

-tp" + (np' + q{x)?p = Xi/j, xe(0,l), (1.6a) 

V'(o) = o, v , (0 = o, (l-6b) 

the potential of which, q{x) = cos<^ s (x), is determined by the known static solution (p s (x). 
Here A = (lu 2 + a 2 /A) 1 ' 2 is the spectral parameter. 

In view of the boundedness of the function < 1 on a finite interval of the 

variable x there exists [H] a countable sequence, bounded from below, of real eigenval- 
ues Amin = A < Ai < A 2 < • • • < A n < . . . of problem (jl.6|) . To each eigenvalue X n 
(nO, 1,2,.. .) there corresponds a single real eigenfunction ip n (x) satisfying the normaliza- 
tion condition 

1 

J i) 2 n {x)dx = 1. (1.7) 


LHere the number of zeroes of the eigenfunction ip n (x) on the interval (0, 1) is equal to the 
index n. In particular, the eigenfunction ipo(x) corresponding to the smallest eigenvalue 
Amin does not have zeroes on (0, 1). 

Since (p(x,p) depends on the set of parameters p, the potential of the SL problem and, 
hence, the corresponding eigenvalues and eigenfunctions of the SL problem also depend 
on the parameters, i.e., A n = X n (p) and tp n = ip n (x,p). 

We shall say [5] that in a certain range of variation of the parameters V C M. the 
static solution (p B (x) is stable in a first approximation with respect to spatiotemporal 
perturbations of the form (jl.5|) if in that region A min (p) > 0. If A min (j9) < 0, then a 
component that is rapidly growing in time appears in the expansion (|1.5|) . and the solution 
(p s (x) is unstable. The points of the vector of parameters which lie on the hypersurface 

Xmin(p) = (1.8) 

in the space V are points of bifurcation (branching) for the solution ip(x,p). The values 
of the parameters for which Eq. (jl.8j) holds are called the bifurcation or critical values for 
the solution (p(x). The corresponding bifurcation curves for the remaining two parameters 
are given by the sections of the surface (jl.8|) by hyperplanes corresponding to fixed values 
of the two parameters. The most important from the standpoint of the possibility of 
experimental verification are the critical curves of the current versus magnetic field type, 

A mi „ (7, h B )=0, (1.9) 
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corresponding to fixed geometric parameters / and a of the junction. 

From a theoretical standpoint knowledge of the bifurcation curves enables one to 
find the number of equilibrium solutions, to understand their structure, and to describe 
the physics of the phenomenon. Numerical simulation simplifies the study and makes it 
possible to estimate the range of variation of the parameters in which one can expect 
stability or instability of the magnetic flux distributions in the Josephson junction. 

Especially important for practical purposes is the ability to check experimentally the 
bifurcation curves of the parameters of the Josephson junction, which in turn is an im- 
portant source of information for refining the physical model. As a concrete example let 
us indicate the methods of studying soliton-like vortex structures of the magnetic flux in 
Josephson junctions on the basis of measurements of the magnetic-field dependence of the 
critical current (see, e.g., Refs. [HI 03 and [7J |HJ HE UD] ) • 

We note that a number of papers have been devoted to the laying of a rigorous math- 
ematical groundwork for reducing the problem of stability of the solutions of nonlinear 
operator equations to an investigation of eigenvalue problems for a linear operator (see, 
e.g., the review ^T] and the collected papers ^2] and the literature cited therein). 

Josephson junctions with an exponentially varying width in the xy plane have been 
studied theoretically and experimentally in Refs. [2] and j3J. The influence of the shape 
on the current-voltage characteristics of the junctions was studied in detail. However, the 
problems involving the determination of the stability regions of the static distributions 
and the structure of the critical curves have not been adequately studied. The present 
paper is devoted to a study of those questions, which are important from the applied and 
theoretical points of view. 

2 Numerical results and discussion 

The exact analytical solutions of equation (jl.3a|) for the case a = are expressed in 
terms of elliptic functions JSI- For o > approximate methods can be used In both 
cases a stability analysis of the solutions using the SL problem ()1.6|) is difficult in view 
of the complexity of the corresponding expressions for the potential q(x). Therefore it is 
advisable to carry out a detailed analysis of the possible states of the magnetic flux in the 
junction and analysis of their stability with the help of a numerical simulation. 

In this paper we use a continuous analog of Newton's method for solution of the 
boundary- value problems ()1.3|) and ()l.b|) . (jl.7|) (see reviews [HJ E3)- The linearized 
boundary-value problems arising at each iteration were solved numerically with the use 
of a spline collocation scheme fTTTj of improved accuracy. 

To permit comparison of our results with the results of Ref. the majority of the 
numerical simulations were carried out for Josephson junctions of length I = 10 and 



For an "infinite" one- dimensional and uniform Josephson junction (<r = 0, I — > oo, 
x G (—00,00)) the well-known exact analytical expression (see, e.g., Refs. |3] and [Hj) 



usually called the fluxon and antifluxon, respectively, is valid. For realistic Josephson 
junctions of finite length those entities are deformed by the geometry of the junction and 
also by the influence of the magnetic field hs and/or the external current 7 and are not 



20. 




(2.1) 
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fluxons in the strict sense of the word (the functions (|2.1|) do not satisfy conditions (jl.2|0 . 
However, because of a number of features of those soliton-like solutions, in particular, 
their finite energy and size, it is appropriate and convenient to use these conventional 
names. For brevity below we denote the fluxon in the Josephson junction as $ 1 . 

According to jH], on the entire axis (—00, 00) the discrete spectrum of the SL problem 
generated by solution (|2.1|) consists of an isolated point A = 0, i.e., the fluxon/antifluxon 
in an "infinite" Josephson junction is found in a quasi-stable (bifurcation) state. 

It follows from general comparison theorems for SL problems that for a finite Josephson 
junction the condition A m i n < holds, i.e., the stability of $ 1 becomes worse. 

For 7 = and small \hs\ the only stable state in a Josephson junction of finite length is 
the Meissner (vacuum) state, which will be denoted by M . For hs = 0, 7 = this "trivial" 
solution of problem (|1.3|) . of the form f(x) = 0, ±27r, ±47r, . . . (there is no magnetic field 
in the junction). For these same parameters the smallest eigenvalue of the SL problem for 
the M distribution is A m ; n (M) = 1. In addition to M there is also an unstable Meissner 
solution (p(x) = ±7r, ±37r, . . ., for which A min = —1. 




(T" 5 10 15 20 5 10 15 20 



Figure 2: Fluxon vortices in a Josephson Figure 3: The ip'(x) curves for different Kb- 
junction. 

For sufficiently large values of the boundary magnetic field hs, stable fluxon vortices 
are generated in the Josephson junction. For example, for Kb — 1.4, 7 = and a = 
0.07, the Josephson junction contains, in addition to the M distribution, the multifluxon 
vortices $ n , n = 1, 2, 3, 4, graphs of which are shown in Fig. |21 The number of vortices is 
determined by the value of the total magnetic flux[5j in the Josephson junction: 

A<p = <p(l)-<p(0). 

For an "infinite" junction (for I — > 00) the value A(p/2ir — > k, where k = 0,2,3, .. . 
is the number of vortices (fluxons) in the distribution (p(x). For the solution in Fig. El 
we have, respectively, Aip($ 1 )/2ir « 1.49, Ay?($ 2 )/27r « 2.48, Ay?($ 3 )/27r w 3.45, and 
A^($ 4 )/27r « 4.36, while for the Meissner solution A^(M)/2vr w 0.49. 

The influence of the external magnetic field Kb on the magnetic flux distribution (p'(x) 
in the junction for the main fluxon $ x at a = 0.07 is demonstrated in Fig. 131 At a certain 
value hs = h m the maximum of the derivative <p'(x) is localized in the middle of the 
junction (curve 2, h B ~ 1.273). For h B < h m the fluxon is "expelled" to the end x = I by 
the "geometric" current g(x) (curve 1, hs = 1)- For hs > h m the fluxon is shifted by the 
external field toward the end x = (curve 2, hs = 1.4). 
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If the length of the Josephson junction is sufficiently large, then a change of the 
current 7, equivalent to a change in magnetic field at the left end x = 0, turns out to have 
only a slight influence on the local maximum magnetic field in the junction, as is well 
demonstrated by Fig. 0] For comparison, the situation in which the shape factor o = is 
demonstrated in Fig. |S] for the same values of the remaining parameters. It is seen that 
variations of the current cause the maximum of the magnetic field to shift to the right or 
left of the center of the Josephson junction, depending on the direction of the current. [TT?] 



1 - 
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Figure 4: The <p'(x) curve for a = 0.07. Figure 5: The ip'(x) curve for a = 0. 

For values Hb > h m the scheme x m of the maximum of the field (f'(x) of the fluxon 
$ x is shifted to the left from the center x = 1/2 toward the x = end (see Fig. EJ). The 
"motion" of the maximum of the function <p'(x) upon variations of the parameters of the 
problem occurs in accordance with the equation 

sin <p(x m , p) + a [<p '(x m , p) - h B ) = , 

which expresses the balance of the Josephson and "geometric" currents at the point x m . 
Here for the $ x distribution the coordinate x m = 1/2 for a = 0, 7 = 0, and any attainable 
hs- This case corresponds to the dashed line in Fig. El which shows graphs of the function 
Xmih-B, cr) for several values of the shape parameter a. Each curve for a > intersects the 
straight line x m (hs, 0) = 1/2 at a certain point h m (a) = tp'(l/2, a) at which the maximum 
magnetic field inside the Josephson junction is centered. It is important to note that for 
cr > the XmiJiB) curves change sharply in the vicinity of h m : slight deviations of the 
external field h b from the value h m cause a significant displacement of the field maximum 
from the center of the junction. 

Let us consider the influence of two geometric parameters of the model, viz., the shape 
parameter a and length /, on the magnetic flux distribution in the Josephson junction. 

The dependence of the smallest eigenvalue of the SL problem for the main fluxon 
on the shape parameter a is shown in Fig. [TJ It is seen that for fixed hs and 7 there exists 
a certain maximum value of a above which the distribution of $ : loses stability, i.e., a 
bifurcation of the vortex occurs upon a change in a. Large values of the magnetic field 
hs correspond to large critical values of a. The value of the current is important at small 
values of a and plays practically no role at values of o close to the maximum. Figure [7| 
thereby demonstrates the destructive character of the shape parameter at large values of 
it and also the stabilizing role of the boundary magnetic field Kb- 
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h B a 

Figure 6: "Movement" of the maximum of Figure 7: Bifurcation upon a change in a. 
the magnetic field in a Josephson junction 
upon a change in Kb- 

The influence of the length I of the Josephson junction on the stability of the main 
fluxon is shown in Fig. |HJ It is seen that A m j n (7) ~ const at I > 12, and therefore to 
a certain accuracy the Josephson junction can be considered "infinite" for $ 1 . At / < 6 
the smallest eigenvalue of the SL problem (jl.fi j) falls rapidly, going to zero at I « 5.23. 
Thus there exists a minimum length of the junction for which the fluxon maintains 
stability. [18 The minimum length clearly depends on all the remaining parameters of the 
model — the external magnetic field hs, the external injection current 7, and also the 
shape parameter a. An analogous statement is also valid for multifluxon distributions 
of the magnetic flux in the junction, including unstable ones. The minimum length for 
multifluxon distributions <3? n falls off rapidly with increasing index n. For example, at 
otherwise equal parameters the vortex <3> 2 exists in Josephson junctions with I > 11.64 
(see Fig. |BJ). Consequently, if a = 0.07, Kb = 1, and 7 = 0, then at lengths / < 5.23 the 
junction is in all respects short, and the only stable distribution in the junction is the 
Meissner one. For 5.23 < I < 11.64 the junction is short for $ 2 , but a nontrivial stable 
vortex $ x can exist in it. 
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Figure 8: Influence of the length of the Josephson junction on thestability of $ x and $ 2 . 

Let us now consider the question of constructing by numerical means the critical 
current versus magnetic field relation, which is determined implicitly by Eq. (jl.iJj) for 
each magnetic flux distribution in the Josephson junction. The importance of this problem 
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stems from the possibility of measuring this relation experimentally. [HJ El El HE EH] We 
note that for different configurations of the magnetic flux in the Josephson junction the 
values of the critical parameters — in particular, the critical current and magnetic field 
- can be substantially different. One should therefore identify the critical parameters for 
specific distributions and for the Josephson junction as a whole. 

The curves of A min (/iB) for the M distribution and the first few stable vortices in a 
Josephson junction at current 7 = and o = are demonstrated in Fig. El For comparison 
the analogous curves for a = 0.07 are shown in Fig. Each curve has two zeroes, the 
distance between which determines the stability region of the vortex upon variation of 
the magnetic field hs- The zeroes themselves are critical values of the field Hb at zero 
current 7. 




Figure 9: The A min (/zs) curves for a = 0. Figure 10: The F{h B ) curves for a = 0.07. 



It is well seen that the role of the shape parameter a is most important at small values 
of Kb- In particular, the $ 1 vortex at a = exists already at h B ~ 0.054. At a = 0.07, 
however, the existence region in respect to field Hb is considerably compressed, and the 
vortex exists starting at h B > 0.75. The amount of compression of the A m i n (/iB) curves 
for different vortices decreases rapidly with increasing hs- 

To complete the picture, Fig. [TT] shows graphs of the total energy of several magnetic 
flux distributions in the Josephson junction 



F(p) 



2<f' 2 + ( x - cos v) 



dx - h B A V - lj(p(0) , (2.2) 



as a function of the external magnetic field hs for a = 0.07 and current 7 = 0. The 
energy values are divided by the total energy Ff^ 1 ] = 8 of an isolated fluxon in an 
infinite Josephson junction. j3] The solid and dashed curves show the energy of the stable 
and unstable distributions, respectively. The points of tangency of the solid and dashed 
curves are the points at which the magnetic flux loses stability. 

Let us now consider the situation for 7 7^ 0. To find the dependence of the critical 
current 7 C on the external magnetic field hs it is necessary to determine the stability 
region with respect to current for the distributions M, $ 1 , $ 2 , etc. The results of such 
calculations for certain values of hs are given in Figs. ^1 — IT41 

Figure H21 shows the A m i n (7) curves for stable solutions of the problem (jl.3|) in a field 
hs = 1.6. The distances between zeroes of the functions are the stability intervals of the 
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Figure 11: The F{h B ) curves for o = 0.07. Figure 12: The A min (7) curves for h B = 1.6. 

corresponding distributions with respect to the current 7. The right and left zeroes of 
the function A m i n (7) are the positive and negative critical currents, respectively, of the 
distribution in the given field hs- Because of the asymmetry of the boundary conditions 
(|1.3b|) and (jl.3cj) for 7 > the critical current of the Meissner distribution [which we 
denote by j c (M)] is the highest, but for 7 < the largest in modulus is the critical 
current 7 C ($ 1 ) of the vortex $ 1 . Consequently, in a field Kb = 1.6 the positive critical 
current of the junction is 7 C = 7 C (M), while the negative critical current is 7 C = 7 C ($ 1 ). 

We note that for the junction geometry under consideration the curve for the M 
distribution has a characteristic plateau — the "breakoff" of the Meissner solution sets 
in at a rather large modulus of the external current. 




y y 



Figure 13: The A m i n (7) curves for $ 1 . Figure 14: The A min (7) curves for <3> 2 . 

Figure [13] demonstrates the A min (7) curves for the main fluxon $ 1 at a = 0.07 and 
several values of Hb- The analogous curves for the $ 2 distribution are shown in Fig. IT41 
We note that with increasing external magnetic field Hb the stability region of the vortices 
with respect to current is narrowed. Consequently, by varying h B one can construct the 
bifurcation curves for individual vortices to acceptable accuracy and from them determine 
the critical values of the current 7 C for a Josephson junction. 

A method of direct calculation of the bifurcation points of the vortices in a Josephson 
junction was proposed in Refs. jTU] and [2Hj . 

Figure shows the critical curves (1.8) for the main stable fluxon $ : for values of 
the shape parameter a = 0, a = 0.001, and a = 0.07. The solid curves correspond to 
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a current 7 > and the dashed curves to 7 < 0. We note that at a value hs ~ 1.273 
the critical curves corresponding to 7 > intersect with each other and with the curve 
corresponding to a = in some narrow region. This means that in that region of hs 
the critical current depends only slightly on the shape of the junction. Geometrically the 
influence of a reduces to a rotation of the critical curves clockwise about the curve for 
ex = by an angle that depends on a . An analogous effect takes place for the critical 
curves corresponding to 7 < 0. 

The critical curve 7 c (/ib) for a junction is constructed as the envelope of the critical 
curves corresponding to different magnetic flux distributions in the junction. In other 
words, the critical curve consists of pieces of the critical curves for individual states with 
the largest module of the critical current at a given Hb- The part of the critical curve 
corresponding to the interval Hb £ [0, 2.8) is illustrated in Fig. For example, let 
Hb = 1-4. At a current 7 = there are five different magnetic field distributions in 
the junction, which are described above (see Fig. |2J). With increasing current 7 in the 
direction of positive values the vortices lose stability in the following order: 

<^> 4 — > <^> 3 — > <^> 2 — ^ <^> x — ^ ikT . 

The last to break off is the Meissner distribution, the critical current of which, j c (M) ~ 
0.156, is the critical current of the junction at a fixed value of the external field. Conse- 
quently, the resistive regime in the junction at hs = 1-4 exists for 7 > 0.156. 

If the current 7 is increased from zero in the negative direction, then the breakoff of 
the distributions occurs in the opposite order: 

M -> 3? 1 -> $ 2 -> $ 3 -> $ 4 , 
and the critical current for the junction will be that of the vortex $ 4 : 7 C ($ 4 ) ~ —0.039. 




Figure 15: The critical curve of the fluxon Figure 16: The critical curve of the junc- 
•D 1 . tion. 

Analogously, for example, at hs — 2.1 the positive critical current for the junction is 
that of the vortex $ 5 , 7 C ($ 5 ) ~ 0.034, while the negative critical current is that of the 
vortex $ 6 , 7 C ($ 6 ) « 0.024, etc. 

Of course, the "switching" of the vortices depends on all the remaining parameters of 
the problem. For example, for a junction with I = 10 and a = 0.07 in an external field 
h B = 1.6 the order of breakoff is as follows: — > $ 2 — > M if 7 > 0, and M — > $ 2 — >• 
if 7 < 0, as is clearly seen in Fig. El 
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We note that our numerically constructed critical curve of the Josephson junction (see 
Fig. Unj) matches well with the theoretical and experimental results presented in Figs. 6a 
and 7a of [3]. 



3 Conclusions 

We have carried out a numerical simulation of the magnetic flux distributions and their 
bifurcations upon variation of the model parameters in a long Josephson junction of 
exponentially varying width. For a stability analysis each distribution is placed in corres- 
pondence with a Sturm-Liouville problem with a potential determined by the given dis- 
tribution. The distributions and, hence, the spectrum of the SL problem depend on the 
model parameters. It is shown that at certain critical (bifurcation) values of one or several 
parameters the magnetic flux in the junction can have a change of stability. For vortex 
configurations in a Josephson junction there are narrow regions of variation of the bound- 
ary magnetic field in which the critical current of the vortices depends only slightly on 
the shape parameter. Corresponding to the magnetic flux vortices are minimum lengths 
of the junction at which the distributions maintain their stability. The critical current 
versus magnetic field curves were constructed by a numerical method for the first few 
stable magnetic field flux distributions in a Josephson junction. The critical curve for the 
junction as a whole is the envelope of the critical curves of the individual distributions. 
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